home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Source Code / Pascal / Applications / NIH Image 1.62b11 / src / Ellipse.p < prev    next >
Text File  |  1996-03-01  |  10KB  |  345 lines
  1. unit Ellipse;
  2.  
  3. interface
  4.  
  5.     uses
  6.         Types, QuickDraw, Palettes, globals, Utilities, graphics;
  7.  
  8.     procedure DrawEllipse;
  9.     procedure GetEllipseParam (var Major, Minor, angle, xxcenter, yycenter: extended);
  10.     procedure ComputeSums (y, width: integer; var MaskLine: LineType);
  11.     procedure ResetSums;
  12.  
  13.  
  14. {Best-fitting ellipse routines by:}
  15.  
  16. {  Bob Rodieck}
  17. {  Department of Ophthalmology, RJ-10}
  18. {  University of Washington, }
  19. {  Seattle, WA, 98195}
  20. {  (206) 543-2449}
  21.  
  22.  
  23. {Notes on best-fitting ellipse:}
  24.  
  25. {  Consider some arbitrarily shaped closed profile, which we wish to}
  26. {  characterize in a quantitative manner by a series of terms, each }
  27. {  term providing a better approximation to the shape of the profile.  }
  28. {  Assume also that we wish to include the orientation of the profile }
  29. {  (i.e. which way is up) in our characterization. }
  30.  
  31. {  One approach is to view the profile as formed by a series harmonic }
  32. {  components, much in the same way that one can decompose a waveform}
  33. {  over a fixed interval into a series of Fourier harmonics over that }
  34. {  interval. From this perspective the first term is the mean radius,}
  35. {  or some related value (i.e. the area).  The second term is the }
  36. {  magnitude and phase of the first harmonic, which is equivalent to the}
  37. {  best-fitting ellipse.  }
  38.  
  39. {  What constitutes the best-fitting ellipse?  First, it should have the}
  40. {  same area.  In statistics, the measure that attempts to characterize some}
  41. {  two-dimensional distribution of data points is the 'ellipse of }
  42. {  concentration' (see Cramer, Mathematical Methods of Statistics, }
  43. {  Princeton Univ. Press, 945, page 283).  This measure equates the second}
  44. {  order central moments of the ellipse to those of the distribution, }
  45. {  and thereby effectively defines both the shape and size of the ellipse. }
  46.  
  47. {  This technique can be applied to a profile by assuming that it constitutes}
  48. {  a uniform distribution of points bounded by the perimeter of the profile.}
  49. {  For most 'blob-like' shapes the area of the ellipse is close to that}
  50. {  of the profile, differing by no more than about 4%. We can then make}
  51. {  a small adjustment to the size of the ellipse, so as to give it the }
  52. {  same area as that of the profile.  This is what is done here, and }
  53. {  therefore this is what we mean by 'best-fitting'. }
  54.  
  55. {  For a real pathologic case, consider a dumbell shape formed by two small}
  56. {  circles separated by a thin line. Changing the distance between the}
  57. {  circles alters the second order moments, and thus the size of the ellipse }
  58. {  of concentration, without altering the area of the profile. }
  59.  
  60.  
  61. implementation
  62.  
  63.     const
  64.         HalfPi = 1.5707963267949;
  65.  
  66.     type
  67.         TMoments= record
  68.                 n: extended;
  69.                 xm, ym,             { mean values }
  70.                 u20, u02, u11: extended; { central moments }
  71.             end;
  72.  
  73.     var
  74.         BitCount, xsum, ysum: LongInt;
  75.         x2sum, y2sum, xysum: extended;
  76.         Moments: TMoments;
  77.         gMajor, gMinor, Theta: extended;
  78.         gxCenter, gyCenter: integer;
  79.         SaveRect: rect;
  80.  
  81.  
  82.  
  83.     procedure DrawEllipse;
  84.  
  85. { basic equations: }
  86.  
  87. {    a: major axis}
  88. {    b: minor axis}
  89. {    t: theta, angle of major axis, clockwise with respect to x axis. }
  90.  
  91. {    g11*x^2 + 2*g12*x*y + g22*y^2 = 1       -- equation of ellipse}
  92.  
  93. {    g11:= ([cos(t)]/a)^2 + ([sin(t)]/b)^2}
  94. {    g12:= (1/a^2 - 1/b^2) * sin(t) * cos(t)}
  95. {    g22:= ([sin(t)]/a)^2 + ([cos(t)]/b)^2}
  96.  
  97. {    solving for x:      x:= k1*y ± sqrt( k2*y^2 + k3 )}
  98.  
  99. {    where:  k1:= -g12/g11}
  100. {            k2:= (g12^2 - g11*g22)/g11^2}
  101. {            k3:= 1/g11}
  102.  
  103. {    ymax or ymin occur when there is a single value for x, that is when:    }
  104. {            k2*y^2 + k3 = 0    }
  105.  
  106.  
  107.         const
  108.             maxY = 1000;
  109.  
  110.         type
  111.             TMinMax = record
  112.                     xmin, xmax: Integer;
  113.                 end;
  114.  
  115.         var
  116.             sint, cost, rmajor2, rminor2, g11, g12, g22, k1, k2, k3: extended;
  117.             xsave, y, ymin, ymax: Integer;
  118.             aMinMax: TMinMax;
  119.             TXList: array[0..maxY] of TMinMax;
  120.  
  121.         procedure GetMinMax (yValue: Integer; var xMinMax: TMinMax);
  122.             var
  123.                 j1, j2, yr: extended;
  124.         begin
  125.             yr := yvalue;
  126.             j2 := sqrt(k2 * sqr(yr) + k3);
  127.             j1 := k1 * yr;
  128.             with xMinMax do begin
  129.                     xmin := round(j1 - j2);
  130.                     xmax := round(j1 + j2);
  131.                 end;
  132.         end;
  133.  
  134.         procedure Plot (x: Integer);
  135.         begin
  136.             MoveTo(gxCenter + xsave, gyCenter + y);
  137.             LineTo(gxCenter + x, gyCenter + y);
  138.             xsave := x;
  139.         end;
  140.  
  141.     begin
  142.         if not EqualRect(info^.RoiRect, SaveRect) then
  143.             exit(DrawEllipse);
  144.         sint := sin(Theta);
  145.         cost := cos(Theta);
  146.         rmajor2 := 1.0 / sqr(gMajor);
  147.         rminor2 := 1.0 / sqr(gMinor);
  148.         g11 := rmajor2 * sqr(cost) + rminor2 * sqr(sint);
  149.         g12 := (rmajor2 - rminor2) * sint * cost;
  150.         g22 := rmajor2 * sqr(sint) + rminor2 * sqr(cost);
  151.         k1 := -g12 / g11;
  152.         k2 := (sqr(g12) - g11 * g22) / sqr(g11);
  153.         k3 := 1.0 / g11;
  154.         ymax := Trunc(sqrt(abs(k3 / k2)));
  155.         if ymax > maxy then
  156.             ymax := maxy;
  157.         ymin := -ymax;
  158.   { Precalculation and use of symmetry speed things up }
  159.         for y := 0 to ymax do begin
  160.                 GetMinMax(y, aMinMax);
  161.                 TXList[y] := aMinMax;
  162.             end;
  163.         xsave := TXList[ymax - 1].xmin;  {  i.e. abs(ymin+1) }
  164.         for y := ymin to ymax - 1 do
  165.             with TXList[abs(y)] do
  166.                 if y < 0 then
  167.                     Plot(xmax)
  168.                 else
  169.                     Plot(-xmin);
  170.         for y := ymax downto ymin + 1 do
  171.             with TXList[abs(y)] do
  172.                 if y < 0 then
  173.                     Plot(xmin)
  174.                 else
  175.                     Plot(-xmax);
  176.     end; { TraceOval }
  177.  
  178.  
  179.     procedure GetMoments;{changed n_}
  180.         var
  181.             x1, y1, x2, y2, xy: extended;
  182.     begin
  183.         with moments, Info^ do begin 
  184.                 if BitCount = 0 then
  185.                     exit(GetMoments);
  186.                 x2sum := x2sum + 0.08333333* BitCount; {NIntegrate[x^2, <x, centerx-0.5, centerx+0.5>]-center^2 = 0.08333333} 
  187.                 y2sum := y2sum + 0.08333333* BitCount; {=correction when the mass of a pixel is seen as an area instead of a point}
  188.                 n := bitcount;
  189.                 x1 := xsum / n;
  190.                 y1 := ysum * PixelAspectRatio/ n;
  191.                 x2 := x2sum / n;
  192.                 y2 := y2sum * sqr(PixelAspectRatio)/ n;
  193.                 xy := xysum * PixelAspectRatio / n;
  194.                 xm := x1;
  195.                 ym := y1;
  196.                 u20 := x2 - sqr(x1);
  197.                 u02 := y2 - sqr(y1);
  198.                 u11 := xy - x1 * y1;
  199.             end;
  200.     end;
  201.  
  202.  
  203.     procedure GetEllipseParam (var Major, Minor, angle, xxcenter, yycenter: extended);{changed n_}
  204. {Return the parameters of an ellipse that has the same second-order }
  205. {moments as those specified by 'm'. }
  206.  
  207. {See Cramer, Mathematical Methods of Statistics, }
  208. {Princeton Univ. Press 1945, page 283.}
  209.  
  210. {The elliptical parameters returned specify an ellipse that}
  211. {has the the same second order moments as that of the profile}
  212. {that generated the moments.  This ellipse need not have the same}
  213. {area as that of the profile, although its area will be close to}
  214. {that of the profile.  In order to refine our measure, we scale}
  215. {the major and minor axes so as to make the area equal to that}
  216. {of the profile. }
  217.   const
  218.    sqrtPi = 1.772453851;
  219.   var
  220.    a11, a12, a22, m4, z, scale, tmp, xoffset, yoffset: extended;
  221.    width, height: integer;
  222.    str1, str2, str3: str255;
  223.    RealAngle: real;
  224.  
  225.  begin
  226.   with info^, info^.RoiRect do
  227.    begin
  228.     width := right - left;
  229.     height := bottom - top;
  230.     if RoiType = RectRoi then
  231.      begin
  232.       major := width / sqrtPi;
  233.       minor := height / sqrtPi * PixelAspectRatio;
  234.       angle := 0.0;
  235.       if major < minor then
  236.        begin
  237.        tmp := major;
  238.        major := minor;
  239.        minor := tmp;
  240.        angle := 90.0;
  241.        end;
  242.       xxcenter := left - 0.5 + width / 2.0;
  243.       yycenter := top - 0.5 + height / 2.0;
  244.       exit(GetEllipseParam);
  245.      end;
  246.    end;
  247.         GetMoments;
  248.         with moments do begin
  249.                 m4 := 4.0 * abs(u02 * u20 - sqr(u11));
  250.                 if m4 <0.000001 then                        
  251.                     m4 := 0.000001;
  252.                 a11 := u02 / m4;
  253.                 a12 := u11 / m4;
  254.                 a22 := u20 / m4;
  255.                 xoffset := xm;
  256.                 yoffset := ym/Info^.PixelAspectRatio;
  257.             end;
  258.         tmp := a11 - a22;
  259.         if tmp = 0.0 then
  260.             tmp := 0.000001;
  261.         theta := 0.5 * arctan(2.0 * a12 / tmp);
  262.         if theta < 0.0 then
  263.             theta := theta + halfpi;
  264.         if a12 > 0.0 then
  265.             theta := theta + halfpi
  266.         else if a12 = 0 then begin
  267.                 if a22 > a11 then begin
  268.                         theta := 0;
  269.                         tmp := a22;
  270.                         a22 := a11;
  271.                         a11 := tmp;
  272.                     end
  273.                 else if a11 <> a22 then
  274.                     theta := halfpi;
  275.             end;
  276.         tmp := sin(theta);
  277.         if tmp = 0.0 then
  278.             tmp := 0.000001;
  279.         z := a12 * cos(theta) / tmp;
  280.         major := sqrt(1.0 / abs(a22 + z));
  281.         minor := sqrt(1.0 / abs(a11 - z));
  282.         scale := sqrt(BitCount * Info^.PixelAspectRatio/ (pi * major * minor ));  {equalize areas }
  283.         major := major * scale;
  284.         minor := minor * scale;
  285.         RealAngle := 180.0 * theta / pi;{force rounding by using real}
  286.         angle := RealAngle;
  287.         if angle = 180.0 then 
  288.             angle := 0.0;
  289.         if major < minor then begin
  290.             tmp := major;
  291.             major := minor;
  292.             minor := tmp;
  293.             end;
  294.         with info^ do begin
  295.                 with RoiRect do begin
  296.                         xxCenter := left + xoffset;
  297.                         yyCenter := top + yoffset;
  298.                     end;
  299.                 SaveRect := RoiRect;
  300.  
  301.             end;
  302.         gxCenter := round(xxCenter);
  303.         gyCenter := round(yyCenter);
  304.         gMajor := major;
  305.         gMinor := minor;
  306.     end;
  307.  
  308.  
  309.     procedure ComputeSums (y, width: integer; var MaskLine: LineType);{changed n_}
  310.         var
  311.             x: longint;
  312.             BitcountOfLine: longint;
  313.             xe, ye: extended;
  314.             xSumOfLine: longint;
  315.     begin
  316.        BitcountOfLine := 0;
  317.        xSumOfLine:= 0;
  318.         for x := 0 to width - 1 do
  319.             if MaskLine[x] = BlackIndex then begin
  320.                   BitcountOfLine := BitcountOfLine+1;
  321.                     xSumOfLine := xSumOfLine + x;
  322.                     x2sum := x2sum + sqr(x);
  323.                 end;
  324.         xsum := xsum + xSumOfLine;
  325.         ysum := ysum + BitcountOfLine * y;
  326.         ye := y;
  327.         xe := xSumOfLine;
  328.         xysum := xysum + xe * ye;
  329.         y2sum := y2sum +  sqr(ye) * BitcountOfLine;
  330.         bitCount := bitCount + BitcountOfLine;
  331.     end;
  332.  
  333.  
  334.     procedure ResetSums;
  335.     begin
  336.         xsum := 0;
  337.         ysum := 0;
  338.         x2sum := 0.0;
  339.         y2sum := 0.0;
  340.         xysum := 0.0;
  341.         bitCount := 0;
  342.     end;
  343.  
  344.  
  345. end.